Analysis of the robustness and dynamics of spin-locking preparations for the detection of oscillatory magnetic fields

Extracting quantitative information of neuronal signals by non-invasive imaging is an outstanding challenge for understanding brain function and pathology. However, state-of-the-art techniques offer low sensitivity to deep electrical sources. Stimulus induced rotary saturation is a recently proposed magnetic resonance imaging sequence that detects oscillatory magnetic fields using a spin-lock preparation. Phantom experiments and simulations proved its efficiency and sensitivity, but the susceptibility of the method to field inhomogeneities is still not well understood. In this study, we simulated and analyzed the dynamic of three spin-lock preparations and their response to field inhomogeneities in the presence of a resonant oscillating field. We show that the composite spin-lock preparation is more robust against field variations within the double resonance effect. In addition, we tested the capability of the chosen composite spin-lock preparation to recover information about the spectral components of a composite signal. This study sets the bases to move one step further towards the clinical application of MR-based neuronal current imaging.


Simulator details
The used rotation matrices ̂( ), ̂( ), ̂( ) can be represented by Each rotation appears counter-clockwise when the axis about which they occur points toward the observer, and the coordinate system is right-handed.
The relaxation matrices , , and can be represented by ∞,̂ is the equilibrium magnetization in the double rotating frame along the spin-lock axis, which we neglect for being much smaller than the generated by the static field 0 . Within the simulator, each rotation step of duration will have two relaxation steps of /2, which means evaluating the four relaxation matrices in = /2. To consider that the target field could oscillate along an arbitrary axis, equation (4) can be generalized as Where and represent the angles formed by the effective field direction with respect to ̂ in the and plane respectively. It follows that when the direction of the target field is along , is 0 and we recover equation (5).
To validate the rotation matrix simulator (RM sim), we compared it with the analytical model presented by Ueda et al. 2018 and the numerical solution of the differential equation presented in this same work using the 4 th -order Rungekutta. The analytical simulator is only valid for a sinusoidal field, so we compared the three simulated signals as a function of the amplitude of the target field. As shown in Supplementary Figure 1. a), both numerical methods approximate well to the Analytical solution. However, the Runge-kutta methods needs a smaller time step (dt < 5E-4 s) to match the accuracy of the RM simulator, as can be seen in panel b). For the same time step dt, the calculation time of the Runge-kutta method is slightly smaller. However, to match the same percentual error, dt must be lower than the required for the RM simulator, making the effective calculation time bigger.
Supplementary figure 1: Validation of the proposed simulation method. a) Signal as a function of the target field amplitude simulated using the analytical solution to the differential equation, the rotation matrix simulator used in this work, and the numerical solution obtained with a 4 th order Runge-kutta method. The Runge-kutta was calculated for different times steps until the error (calculated as the mean distance to the analytical solution) was lower than 1%.

Dynamics in the double rotating frame
The double rotating frame is usually employed to depict magnetization dynamics under SL preparation techniques. The double rotating frame ( ', ', ') is defined as if the simple rotating frame is now also rotating around the SL axis at the SL frequency. We can decompose the target oscillating field ( ) = ( + )̂ into two fields of half the amplitude oscillating in the plane with opposite frequencies as shown in Supplementary Figure 2 i): When in resonance with the SL frequency, one of the fields will be static in the double rotating frame while the other one will be oscillating at 2 . Applying the rotating wave approximation (firstly proposed in this context by Ueda et al., Journal of Magnetic Resonance, 2018), we discard the effect that the off-resonance field will have on the magnetization. Supplementary Figure 2 ii), iii) and vi) shows the magnetization dynamics for the three presented SL preparations in the double rotating frame. For BASL there is only one SL frequency, therefore the target field ( ) remains static and deviates the magnetization from the SL axis. For RESL and CRESL, the same happens during the first half of the SL pulse. When the SL changes sign, the opposite field ( ) will be on resonance, and its initial position will depend on the phase accumulated during the first part of the evolution. Depending on this phase, the contrast accumulated during the first SL part can be compensated. iii) and iv) BASL, RESL and CRESL magnetization dynamics represented in the double rotating frame respectively. As the double rotating frame rotates around the SL axis, the sign of rotation changes when SL1 changes to SL2 in the RESL and CRESL cases. Therefore, only 1 and 1 are present during the first SL part, and 2 and 2 are present during the second part for RESL and CRESL. The magnetization path follows the same color code as the one used in Figure 1 of the main text. The 180° pulse dynamic is indicated with dashed lines.

SL pulse calibration
An error in the excitation angle can be corrected by the change in the sign of the second half of the spin-lock pulse (RESL) and an error in the 0 field can be corrected by the 180° pulse (CRESL). However, no compensation can be made for the 1 inhomogeneity influence during the application of the SL pulse. This will induce a different from the one of the target field and the resonance condition will not be fulfilled. To solve this problem, an SL calibration can be made to find the 1 error for each target frequency. Supplementary Figure 3 a) shows the measurements at different resonant frequencies as a function of and b) the obtained linear regression. The initial phase was set to 0, and 30 repetitions were acquired and averaged for each point on each curve. From the fitted curve, the SL amplitude was corrected before every measurement.

Supplementary figure 3. SL amplitude calibration. a) Percentual change of the signal as a function of for different values of the target oscillating field frequency ( ). b) Linear regression of the detected ( ) vs the expected frequency (
). The obtained slope was then used to define the amplitude of each SL pulse.

BASL contrast in the presence of field inhomogeneities
To explain the result showed in Figure 3 (b.i) of the main text, we fitted the BASL contrast using the field inhomogeneity as fitting parameter. The excitation angle , the off-resonance frequency Δ 0 and the spin locking frequency were used as fitting parameters. The rest of the parameters were set as stated in the methods. Supplementary Figure 4 shows the fitting result. The fitted parameters were = 73°, Δ 0 = 0 Hz, and = 90,88 Hz. This indicates that 1 is the main contributor to the behavior observed in this case. This behavior comes from the magnetization without a target field ( 0 ), when the tip down pulse angle is smaller than 90°. This generates an oscillation around the SL axis that result on a transversal magnetization after the tip up pulse.

SAR estimation
The specific absorption rate (SAR) limits the range of detectable frequencies for the SIRS approach, especially at high field strengths. Previous reports have calculated the SAR of the CRESL sequence; however, the values are highly dependent of sequence definition, coil used and subject scanned. Therefore, to assess the detection limits for our own experimental set-up, we used the SAR check computed by the scanner. In our case, we used a 64-channel head birdcage coil and set the subject values to typical values of high (170 cm) and weight (60 kg). The minimum is limited by . Using the maximum allowed value, = 100 ms, we would have the highest temporal resolution and worst SAR case for the minimum possible = 165 ms (minimum possible value given by preparation duration plus duration of the EPI readout). For this value, the obtained SAR as a function of is shown in Supplementary Figure 5. The SAR value is the same for BASL and RESL, only differing with CRESL by the 180° pulse. The greatest contribution to SAR is given by the SL pulses, therefore it does not differ so much between the proposed preparations. In this setup, the frequency limit before exceeding the limit set by the International Electrotechnical Commission (IEC) (3.2 W/Kg for the head) is 222 Hz, 220 Hz and 310 Hz for BASL and CRESL and the (CRESL)/ approach respectively. If a higher frequency value must be used as target, needs to be increased with the consequent loss of temporal resolution.

Measurement of individual signals
Supplementary Figure 6 shows the point-to-point division between the (CRESL) and acquisition of the 4 composite signals presented in Figure 7 a) of the main text. The signal shows traces of the 0.1 Hz oscillation, despite the acquisition being performed for equal to 30, 60 and 90 Hz. This influence is greater at 30 Hz and diminishes when increasing . This shows that slow drifts cannot be fully filtered out using the ppd between on and off acquisitions. The reason behind this is that the two acquisitions have a time difference of . Therefore, any signal variation that occurs within a will be registered differently in the and signal. An alternative to this is applying a high-pass filter to eliminate the influence of slow drifts if there is prior knowledge of the signals present. Nevertheless, the 90 Hz signal only has 10% of the amplitude of the slow-wave and is still the biggest detected component at = 90 Hz.